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SPATIAL MORAN MODELS 

I. STOCHASTIC TUNNELING IN THE NEUTRAL CASE 

By Richard Durrett 1 and Stephen Moseley 2 
Duke University 

We consider a multistage cancer model in which cells are arranged 
in a d-dimensional integer lattice. Starting with all wild-type cells, 
we prove results about the distribution of the first time when two 
neutral mutations have accumulated in some cell in dimensions d > 2, 
extending work done by Komarova [Genetics 166 (2004) 1571 1579] 
for d = 1. 


1. Introduction. The accumulation of mutations is important not only 
for cancer initiation, progression, and metastasis, but also for the emergence 
of acquired resistance against chemotherapeutics, radiation therapy, or tar¬ 
geted drugs. For this reason there is a large and growing literature on the 
waiting time Tfc until some cell has acquired k prespecified mutations. In all 
the models we consider, type i individuals mutate to type (i + 1) at rate 
Ui- |_i. The dynamics considered have most often been studied in multi-type 
Moran models with a homogeneously mixing population of constant size. 
Here we will concentrate on how results change when one considers a spatial 
Moran models, and as is the case for much earlier work we will concentrate 
on the behavior of T 2 . 

We suppose that cells of type 0 and type 1 have relative fitness 1 and A. 
Since we will only consider the waiting time for the first type 2, the relative 
fitness of type 2’s is not important. In this work we will consider situation 
in which A is so close to 1 that the mutations are essentially neutral. For 
cancer applications, this is a restrictive assumption, and it will be removed 
in the companion paper (part II) by Durrett, Foo and Leder [6]. However, 
the current result applies to the important case of tumor suppressor genes. 
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In that case, when both copies of the gene are inactivated trouble develops, 
but while there is one working copy the cell can function normally. 

We begin by recalling results for the Moran model in a homogeneously 
mixing population of size N. Here and in what follows the mutation rates rq 
and selection coefficient A depend on N, even though this is not indicated 
in the notation, and we write a at <C 6tv if ajv/6 jv —> 0 as N —> oo. The next 
result made its first appearance on page 16,230 of Nowak et al. [17]. Since 
then it has appeared in print a number of times: [10, 11, 14, 18], and in 
Nowak’s excellent book [16] on Evolutionary Dynamics. 


Theorem 1. In the neutral case of the Moran model, A 
sume that 


( 1 ) 


—iv< — 

y/U2 Ul 


and let U\,U 2 —>• 0 then we have 

P(t 2 > t/Nu\uy 2 ) —> exp(— t). 
The same conclusion holds if |A — 1| <C mj 2 . 


1, if we as- 


Durrett and Schmidt [7] applied these ideas to study regulatory sequence 
evolution and to expose flaws in Michael Behe’s arguments for intelligent 
design. Durrett, Schmidt and Schweinsberg [8], see also Schweinsberg [19], 
generalized this result to cover Tfe. 

The conditions in the result may look mysterious but they can be derived 
by simple reasoning. Here and throughout the paper and f(u) ~ g(u) means 
f(u)/g(u ) —>• 1 as u —>• 0. Suppose first that A = 1. 

(Al) If we start the Moran model with fc<JV type l’s and the rest type 
0, then the l’s behave like a critical branching process. The time needed for 
the l’s to die out is 0(k ) and the number of type-1 births before they die 
out is 0{k 2 ). Thus we expect the first type 2 to occur in a type-1 family 
that reaches size k\ = 0(1/ f/uf), and hence has 0(k 2 ) = 0( 1 / 7 x 2 ) births. 
The condition 1 /y/ui <C N in (1) guarantees k\ <C N. 

(A2) Since the voter model is a martingale, the probability a type-1 mu¬ 
tation creates a family that reaches size 1 /\/u 2 is yjui- More to the point a 
simple computation (consider what happens at the first jump) shows that 
the probability a type-1 family gives rise to a type 2 before it dies out is 
~ y/u 2 ■ Since mutations to type 1 occur at times of a rate Nui Poisson pro¬ 
cess and with probability ~ give rise to a type 2, it follows that if p 2 is 
the birth time of the type-1 mutant that first gives rise to a type 2 then 

P(p2 > t/Nu\uy 2 ) -A exp(— t). 

To complete the proof we need to show that T 2 — P 2 ^ P 2 , and for this we 
need the condition N <C 1/ui in (1). 
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(A3) By the discussion of (Al), the first mutation will occur in a family 
that reaches a size 0(1/y/u^). If [A — 1| <C u 2 , then computations with 
Girsanov’s formula show that (in the limit U 2 —>• 0) the behavior of the 
Moran model, while it is 0(1 is indistinguishable from the case with 
no drift. 

The assumption of a homogeneously mixing cell populations simplifies cal¬ 
culations considerably, but is not realistic for solid tumors. For this reason, 
Komarova [12] considered a spatial model, which is very similar to one intro¬ 
duced much earlier by Williams and Bjerknes [20]. Due to work of Bramson 
and Griffeath [2, 3], the second model is known to probabilists as the biased 
voter model. 

In the usual formulation of the biased voter model, each site on the 
d- dimensional integer lattice 7L d can be in state 0 or 1 indicating the presence 
of a cell with relative fitness 1 or A > 1. Cells give birth at a rate equal to 
their fitness, and the offspring replaces a nearest neighbor chosen at random. 
When A = 1 this is the voter model which was introduced independently by 
Clifford and Sudbury [4] and Holley and Liggett [9] . For a summary of what 
is known see Liggett [15]. 

In the biased voter model births drive the process. In Komarova’s version 
cells die at rate 1 and are replaced by a copy of a nearest neighbor chosen 
with probability proportional to its fitness. A site with m neighbors in state 
i makes 

transitions at rate 

0 —y 1 Ani/(Ani + no) 

1—^0 n 0 /(Ani + n 0 ) 

In d = 1 if the set of sites in state 1 is an interval [l, r] with i < r then any 
site that can change has ni = no = 1 so Komarova’s model is a time change 
of the biased voter model. In d > 2 this is not exactly true. However, we are 
interested in values of A = 1 + s where s = 0.02 or even less, so we expect 
the two models to have very similar behavior. In any case, the difference 
between the two models is much less than their difference from reality, so 
we will choose to study the biased voter, whose duality with branching 
coalescing random walk (to be described below) gives us a powerful tool for 
doing computations. 

Since we want a finite cell population we will restrict our process to be a 
subset of (—L/2,L/2] d . Komarova [12] uses “Dirichlet boundary conditions”, 
that is, she assumes her space is an interval with no cells outside, but this 
is awkward because the set of type-1 cells may reach one end of the interval 
and then no further changes happen at that end. To avoid this, we will use 
periodic boundary conditions, that is, we consider (Z mod L) d . The resulting 
toroidal geometry is a little strange for studying cancer. However, using 
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(Z mod L) d has the advantage that the space looks the same seen from any 
point. Our results will show that for the parameter values the first type 2 
will arise when the radius of the set of sites occupied by l’s is <C L so the 
boundary conditions do not matter. 

Let be the set of cells equal to 1 in the voter model with no mutations 
from 0 to 1 on 7L d starting from a single type 1 at 0. Let |£^| be the number 
of cells in and let 

(2) v d = 1 -Eexp(^-u 2 J |£° | ds 

This quantity, which is defined for the voter model without mutation, cal¬ 
culates the probability, which depends on the dimension d, that a mutation 
to type 1 gives rise to a type 2 before its family dies out. To see why this is 
true note that the integral / Q T ° |£°|cZs gives the total number of man-hours 
in the type -1 family, and conditional on this the number of mutations that 
will occur is Poisson with mean u 2 f |^| ds. 

Since mutations to type 1 in a population of N cells occur at rate Nu\ 
this suggests that 


( 3 ) 


P{t 2 >t)-> ex.p(—Nuii>dt). 


As we will explain in a moment, there is a constant 7 d so that v d r^j ldhd{u2) 
as U 2 —>• 0 where 


( 4 ) 


{ u 1 / 3 , 

u 1 / 2 log 1 / 2 (l/u), 
u l /\ 


d= 1, 
d = 2, 
d> 3. 


To state the result we need one more definition: 

tt 1 / 3 , d = 1, 

(5) gd(u) = ^ log - 1 / 2 (l/u), d = 2, 

1, d> 3. 


Theorem 2. In the neutral case of the biased voter model, A = 1. if we 
assume 

1 «jv«»A2) 


/g\ 

h d {u 2 ) ui 

then there are constants 7 ^ given in (12) and (13) so that as ui,U2->-0 
P(t 2 > t/Nui'y d h d (u 2 )) —t exp(-t). 

The same conclusion holds if |A — 1| <C h d [u 2 ). 
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In d = 1 this result was proved by Komarova [12], see her equation (62) 
and assumption (60), then change notation u\ -A u, U 2 -A u\. See also her 
survey paper [13]. Note that when d > 3 the order of magnitude of the 
waiting time and the assumptions are the same as in Theorem 1. In d = 2 
there are logarithmic corrections to the behavior in Theorem 1, so only in 
the case of d = 1 (which is relevant to cancer in the mammary ducts) does 
space make a substantial change in the waiting time. 

The reasons for the conditions in Theorem 2 are the same as in Theo¬ 
rem 1. 

(Bl) We will see that the mutation to type 2 will occur in a type-1 family 
that reaches size k = 0(l/hd(u2 ))• The left-hand side assumption in (6) 
implies that k -C N, so the type-2 mutant arises before the l’s reach fixation. 

(B2) Let p 2 be the time of the first type-1 mutation that begins the family 
that eventually leads to a type 2. Since mutations to type 1 occur at rate 
Nu\ and lead to a type 2 with probability I'd, it is easy to see that 

P(p 2 >t)—> exp(-Nuii'dt) 

so to prove the result we need to show that with high probability 72 — 
P 2 ^ P 2 - As the reader will see, this is guaranteed by the right-hand side 
assumption in (6). 

(B3) As in the discussion of Theorem 1, once we know that the mutation 
to type 2 will occur in a type-1 family that reaches size k = 0 (l/hd(u 2 )), it 
follows that if [A — 1| <C hd{u 2 ) then (in the limit U 2 -A 0) the behavior of 
the size of the biased voter |£°| is, while it is 0 (l/hd(u 2 )), indistinguishable 
from the case with no drift. 


2. The key to the proof. The size of the voter model, when |£j? | > 0, is 
a time change of symmetric simple random walk, with jumps happening at 
two times the size of the boundary |<9£°|, which is the number of nearest 
neighbor pairs with and y ^ The one-dimensional case is easy 

because when $ ^ 0 the boundary |dtf \ = 2. The key to the study of the 
process in d > 2 is the observation that there are constants Pd so that 


(7) 


m, 


(2dp d \$\, d> 3, 

l4/3 2 |e t °i/iog(ie t °i), d= 2 , 


where |<9£°| /(|£?|) means that when [£°[ is large, |//(|£t ) |) is close to 

1 with high probability. 

The intuition behind this result is that the voter model is dual to a collec¬ 
tion of coalescing random walks, so in d > 3 neighbors of points in will be 
unoccupied with probability ~ Pd, the probability two simple random walks 
started at 0 and ei = (1,0,... ,0) never hit. In dimension d = 2, the recur¬ 
rence of random walks implies that when |£°| = k is large, most neighbors of 
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points in will be occupied, but due to the fat tail of the recurrence time 
sites will be vacant with probability ~ fa /log k, where fa = it. 

Before we try to explain why (7) is true, we will list an important conse¬ 
quence. Let Tfc be the first time \f/\ = k. Let 

( n 2 , d= 1, 

a n = < 2nlogn, d = 2, 

[n, d> 3. 


Lemma 1. Let be the unbiased voter model (i.e., \ = 1) starting from 
a single occupied site. 


( 8 ) 


ie° 


'I'ne’^'CLnt' 


n 


T ne < oo => (Y t \Yo = e), 


where => indicates convergence in distribution of the stochastic processes and 
the limit has 


dY t = 


V2dB t , 

VWdYtdBt, 


d = 1. 

d> 2, 


where Bf is a one-dimensional Brownian motion. In d = 1 the process is 
stopped when it hits 0. In d>2, 0 is an absorbing boundary so we do not 
need to stop the process. 


In d = 1 this result is trivial. If one accepts (7) then (8) can be proved 
easily by computing infinitesimal means and variances and using standard 
weak convergence results. In d> 2, (7) and (8) are almost consequences of 
work of Cox, Durrett and Perkins [5]. They speed up time at rate a n , scale 
space by 1 / y/dff and assign each point occupied in the voter model mass 
1/n to define a measure-valued diffusion X n which they prove converges 
to super-Brownian motion. See their Theorem 1.2. (Their scaling is a little 
different in d = 2 but this makes no difference to the limit.) 

Let V/ S {x) be the fraction of sites adjacent to x in state 0 at time s (with 
the prime indicating that we multiply this by logn in d = 2, see page 196). 
A key step in the proof in [5] is to show, see (II) on page 202, that for nice 
test functions <f> 


(9) 



—> 0 , 


where X/(f ) denote the integral of the function / against the measure Xf. 
The result in (9) shows that when we integrate in time and average in space 
(multiplying by a test function to localize the average) then (7) is true. 

From the convergence of the measure valued diffusion X n to super- 
Brownian motion, (8) follows by considering the total mass. Earlier we said 







TUNNELING IN SPATIAL MORAN MODELS 


7 


(8) is almost a consequence of [5], since they start their process from an ini¬ 
tial measure [i.e., O(n) initial l’s] while consider a single occupied site and 
condition on reaching ne. However, this defect can be remedied by citing 
the work of Bramson, Cox and LeGall [1], who have a result, Theorem 4 on 
page 1012 that implies (8) in d> 2. 

The result in (8) is enough for Section 3, but for the calculations in Sec¬ 
tion 4 we will need a version of (7). In that section we will compute under 
the assumption that if |£f| = k 


( 10 ) 


incOi _\ 2d Pdk, d> 3, 

1 ^ l ~\^ 2 k/logk, d = 2. 


If one wants to give a rigorous proof of the estimates there, then small values 
of k, can be treated with the inequalities 

Ck 1/d <\d$\<2dk, 


and one can control large values of k using (9) and estimates such as (Jl) 
and (J2) on page 208 of [5]. We will assume (10) in order to avoid getting 
bogged down in technicalities. 


3. Proof, part I. Let v e d be the probability defined in (2) ignoring mu¬ 
tations to type 2 that occur before T ne . The size of the voter model, |£°|, is 
a martingale, so if we let P\ to denote the law of the voter model starting 
from one occupied site P\(T n£ < oo) = 1/ne. Applying (8) now, 


( 11 ) 




ne 


rTo 

1 — E e exp ( —na n U 2 / Y s ds 


where To = min{f: Y t = 0}, E e is the expected value for (Y t \Yo = e). We have 

f n 3 , d = 1, 

na n =< 2n 2 logn, d = 2, 

l n 2 , d> 3. 


So if we set n 


l/h d (u 2 ) then (4) implies na n U 2 —» 1 and using (11) gives 


v e d ~ h d (u 2 ) 


' 1 - E e exp(— J 0 To Y s ds)' 
£ 


Thus the type-2 mutation will occur in a family that reaches sizes 0(l/h d (u2)), 
and we must assume l/h d (u 2 ) N. 

If we ignore the time to reach size l/h d (u 2 ), the time needed to generate 
the type-2 mutation is, by (8), of order 


a(l/h d {u 2 )) ~ 


u, 


-2/3 


2n 2 1/2 log 1/2 (l/n 2 ), 


u. 


- 1/2 


d = 1, 
d = 2, 
d > 3, 
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where we have written a(n) for a n for readability. Thus for (B2) we need 
a(l/hd(u 2 )) <C l/N Ul h d (u 2 ), which means N <C 9 d(u 2 )/ui. 

The next order of business is to compute v&. Stochastic calculus (or cal¬ 
culations with infinitesimal generators) tells us that 

/ rT 0 

v(x) = E x exp I — / Y s ds 

is the unique function on [0, oo) with values in [0,1], u(0) = 1 and 
v" — xv = 0 in d=l. (3dxv" — xv = 0 ind>2. 

In d = 1 all solutions have the form: 


v(x) = aAi(x ) + (3Bi(x), 


where Ai and Bi are Airy functions 



dt. 


Since Bi is unbounded and Ai is decreasing on [0,oo), we take ft = 0 and set 
a = 3 2 / 3 T(2/3) to satisfy the boundary condition, u(0) = 1. Letting e —>• 0 
we conclude that 


(12) 7i = -aAi'{ 0) = 3 1/3 T(2/3)/T(l/3). 

_ \l‘2 l 

In d > 2, v{x ) = exp(— /3 d ; x), and we have 

(13) 7 d = Pd 1/2 - 


4. Proof, part II: Missing details for A = 1. In the previous section we 
have calculated the probability v e d that a type-1 family reaches size e/hd(u 2 ) 
and then gives rise to a type 2. To let e —> 0 and prove Theorem 2 we need 
to consider the possibility of a mutation to type 2 in a family that (i) never 
reaches size ne, or (ii) will reach ne but has not yet. To have a convenient 
name we will call these small families. Families of the first kind arise at rate 
Nu\(l — 1/ne) and families of the second kind arise at rate Nui/ne. We 
will now calculate the expected rate at which type 2’s are born from these 
small families. In the proof of Theorem 2, we will let £ —>• 0 slowly as n —>• oo 
so we can and will assume ne —> oo. 

Consider the voter model starting from a single 1 at the origin at time 
0. Let 14 be the total time spent at level k, that is, \{t : \^\ = k}\ and let 
be the total number of returns to level k before leaving the interval (0,ne). 
Recalling our assumption in (10), we let q{k) the rate jumps occur at level k. 
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Let Sk be the embedded discrete time chain, which is a simple random 
walk, and let T7 1 " = min{n > 1: S n = k}. 


Ei (^ T ° |€°| ds T 0 < T n ^j = E\ ^ kV k 


Ti < T n 


(14) 


My 


kN k 


To < T ne 


.ti 

vA Ti(Tfc<oo) k 
^P fe (T+>T 0 )g(fc) : 


where the bar indicates conditioning on To < T ne . A similar argument shows 
that 

/•T \ ne: 

(15) 


/ r-Lne \ 

Hi 


k 


> T ne ) q(k)' 

where the hat indicates conditioning on T n£ < To. 

The three conditional probabilities we need can be computed using facts 
about simple random walk that follow from the fact that it is a martingale. 

n frr „ n Pi(T k < oo)P k (T 0 < T ne ) (1/k )(1 - k/ne) 

(16) Piffi,<oo) =- Pl(To<r „ e) -= (i — 1/ne) ' 

For the next two we note that the first step has to be in the correct direction 
for these events to happen. 

( 1 / 2 )( 1 /k) 


(17) 


(18) 


P k (T+>T 0 ) = 
T fc (T+ > T ns ) = 


(1 — k/ne) ’ 

(1/2) (1 /(ne — k)) 
(k/ne) 


Thus the expected total man-hours f^° | ds for a family that will die out 

before reaching size ne is 


(19) 


2 Ti h 

£(1 - k/ne)° 


(1 — 1/ne) 


k =1 


q(k)' 


and in families that have yet to reach size ne, 

, , 2 , k 2 

(20) — ^(ne-fc) 


fc=i 


q(k)' 


The next result shows that the contribution of small families are indeed 
negligible. Note that in all three cases the order of magnitude of the contri¬ 
butions from small families is the same as the overall rate, but contains a 
constant that —>• 0 as e -A 0. 
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Lemma 2. The expected total man-hours in small families is 


< 


N Ul u 2 /3 • 


Nmu 2 2 log 1//2 (l/u 2 ) • 


7e 

24^’ 


Nu\u 2 


1/2 


2 dp d ' 


d= 1, 
d = 2, 
d> 3. 


Proof. In one dimension, q(k) = 2. The sum in (19) is dominated by 

^ r n £ (nr-\2 


f 


(1 — xjneyxdx = 


(ne) 2 


rn£ 2 , , i (ne)- 

y (ne - y)dy = —py 


to 


Thus, families of the first kind produce type 2’s at rate < -/VuiU2(ne) 2 /12. 
The expression in (20) is dominated by 


2 f n£ 

— / (ne — x) 
ne Jo 


-x)Jdx= {ncf 


6 


Thus, families of the second kind produce type 2’s at rate < Nu\U 2 (ne) 2 /6. 
Adding the last two conclusions gives the result for d = 1. 

In d > 3, (10) implies q(k) = 2 d/3dk, so (19) becomes 


1 

dPd 


£(1 - k/nef. 


k =l 


The sum is bounded above by the integral 

rns 

/ (1 — x/ne) 2 dx= , 

Jo 3 


ne 


so with our choice of n = u 


- V 2 


families of the first kind produce type 


1 /2 

2’s at rate bounded above by Nu\u 2 e/(3d/3d)- Setting q(k ) = 2 dfidk, (20) 
becomes 


^(ne — k)k. 


df3dne 


k =l 


The sum is bounded above by the integral 


f 


(ne — x)xdx = 


(ne) 3 


1 /2 

Thus, families of the second kind produce type 2’s at rate < Nu\u 2 e/(6 d/3d)- 
Adding the last two conclusions gives the result for d> 3. 
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In d = 2, (10) implies q(k) = 4/?2 k/logk, so (19) becomes 


2/3 2 


Ed — k/ne) 2 log k. 


k =1 


Each term in the sum is bounded above by log(ne), so the sum is less than 
nelogne. Since n = u 2 1/7 “ \og~ 1 ^ 2 (l/u 2 )-, families of the first kind produce 
type 2’s at rate bounded above by 

Nuiu 2 • -^-nelog(ne) = Nuiu 2 ■ 1/2 log -1/2 (l /u 2 ) ■ - log(l/tt 2 ) 

Zp 2 Zp 2 Z 

= j^Nu 1 uY 2 \og 1,2 {l/u 2 ). 

Taking q(k) = Afok/logk, (20) becomes 

\ ( ne \ 

^(ne — k)k\ogk 


2 /32ns \ 

The sum is bounded above by 


,k =1 




f 

■JO 


(ne — x)xlog(ne) dx < - -log(ne). 

6 


Thus families of the second kind produce type 2’s at rate bounded above by 
Nu\u 2 1 


ne 2/32 ne 6 


(ne')'^ 1 

log(ne) = Nu\U 2 ■ nelog(ne) 


12/3 2 


24/3 2 


Nuiu^ 2 log 1 / 2 (l/« 2 ). 


Adding the last two conclusions gives the result for d = 2 and completes the 
proof. □ 


5. Proof, part III: Almost neutral mutations. In the biased voter model, 
whose law we denote by P x , jumps occur at rate 1 + A times the size of the 
boundary. To compensate for this we need to run the unbiased (A = 1) voter 
at rate (1 + A)/2. If we do this, call the resulting law P°, and let lot is a 
realization of $ run up to time T then the Radon-Nikodym derivative 


dP x 

dP° 


(lot) 




where n+ and n_ are the number of up jumps in oot when 0 < t < T. 

If maxt^T \£t\ = O(K) then the difference 0 < n + — n_ = 0(K). Since 
under Pq, |£f| is a time change of simple random walk, we see that the total 
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number of jumps n+ + n_ = 0(I\ 2 ). Taking K = 1//id(if 2 ) and assuming 
|A — 1| <C hd(u 2 ), when U 2 is small the Radon-Nikodym derivative is 


A — 1 \ n+ 
+ A + l) 


1 - 


A — 1 


= 1 + 


A- 1 
A + l 


n_|_—n_ 


A + l, 

(A - l) 2x n++n ~ 


1 - 


(A + l) s 


The last result implies that (8) extends to almost neutral mutations, and 
the computations in Section 2 are valid. To extend the part of the proof 
in Section 3, we need to check that (16)-(18) are true asymptotically for 
almost neutral mutations. To do this we recall that if a < x < b 

(21) P x \T b <T a ) = 9 ^ where 6 = 1/A. 


When 0<a<x<6 = 0{l/hd{u2)) and |A — 1| <C hd(u 2 ) we have 


P*{T b <T a 


x — a 
b — a 


To show that the sums come out the same we need the following uniform 
version which follows from (21). If |A — l\hd(u 2 ) —>• 0 then for any C fixed 

P*(T b <T a ) 


sup 

0<—a,b<C/h c i(u2) 


—a/ (b — a) 


- 1 


0 . 
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